clear all 
set more off 
set maxvar 15000 
clear matrix

/* Bring in multiple imputation output, 
   which includes baseline */
 
    use "$Mydirectory2/appendix_b/MI_output_post1940.dta", clear
    replace decade=decade-2 if estimate==2
    
    keep est_* estimate decade 
    
    tempfile output1
    save `output1'
    

* Bring in output from full distribution 
    use "$Mydirectory2/appendix_b/beta_IGE_output.dta", clear
    
    gen est_=. 
    gen est_ub_=. 
    gen est_lb_=. 
    
    //Save means and confidence intervals
    forval d=1/7 { 
        sort iter 
        
        * Grab mean
        sum beta_`d'
        replace est_=`r(mean)' if _n==`d' 

        * Now grab confidence intervals using percentiles
        xtile beta_perc = beta_`d', nq(100)
            
        sum beta_`d' if beta_perc>=97 & beta_perc<=98   
        replace est_ub_=`r(mean)' if _n==`d' 

        sum beta_`d' if beta_perc>=2 & beta_perc<=3 
        replace est_lb_=`r(mean)' if _n==`d' 
        
        drop beta_perc
    } 

    /*Save relevant variables in right 
      format for append below */
    keep est_*
    keep if est_!=.
    
    gen estimate=3
    gen decade = 1900+(_n*10)
    
    tempfile tempie
    save `tempie'
    
* Append 2 files
    use `output1', clear
    append using `tempie'
    
    sort decade estimate
    
    replace decade=decade+1.5 if estimate==2
    replace decade=decade+3 if estimate==3
    
    #delimit ;
    twoway (scatter est_ decade if estimate==1,  msymbol(circle) mcolor(midblue) msize(medium) yaxis(1)) 
           (rcap est_lb_  est_ub_  decade if estimate==1, lpatter(solid) lcolor(midblue) yaxis(1) lwidth(0.5))
           (scatter est_ decade if estimate==2,  msymbol(triangle) mcolor(green*0.75) msize(small) yaxis(2)) 
           (rcap est_lb_  est_ub_  decade if estimate==2, lpatter(solid) lcolor(green*0.75) yaxis(2) lwidth(0.5) )
           (scatter est_ decade if estimate==3,  msymbol(square) mcolor(purple*0.75) msize(small) yaxis(2)) 
           (rcap est_lb_  est_ub_  decade if estimate==3, lpatter(solid) lcolor(purple*0.75) yaxis(2) lwidth(0.5) )
           ,
    xti(" " "Decade of respondent's birth") xlabel(1910(10)1970) xscale(range(1905 1975))
    ylabel(0(.25)1, axis(1)) yti("IGE estimate (baseline)" " ", axis(1)) scheme(s1color) 
    ylabel(0(.1)0.3, axis(2))  yti(" " "IGE estimate (alternative imputations)", axis(2)) 
    legend(on ring(0) size(medsmall) row(3) pos(8) order(1 "Baseline" 3 "Multiple imputations" 5 "Simulation, full dist." ))
    xlabel(1910 "1910s" 1920 "1920s" 1930 "1930s" 1940 "1940s" 1950 "1950s" 1960 "1960s" 1970 "1970s", labsize(small) ) ;  
    #delimit cr
    graph export "$Mydirectory2/appendix_b/MI_results.pdf", as(pdf) replace

    